 ///////////////////////////////////////////////////////////////////////////////
//! @file	: integral.cl
//! @date   : 2016/04/26
//! @author: guyadong
//! @brief  : Calculates the integral matrix of an image
////////////////////////////////////////////////////////////////////////////////
#include "common_types.h"
#ifndef CL_DEVICE_LOCAL_MEM_SIZE
#error not defined CL_DEVICE_LOCAL_MEM_SIZE by complier with options -D
#endif
#ifndef SRC_TYPE 
#error not defined SRC_TYPE by complier with options -D
#endif
#ifndef DST_TYPE 
#error not defined DST_TYPE by complier with options -D
#endif
#ifndef INTEGRAL_COLUMN_STEP
#error not defined INTEGRAL_COLUMN_STEP , must be 2,4,8,16
#elif (INTEGRAL_COLUMN_STEP!=16)&&(INTEGRAL_COLUMN_STEP!=8)&&(INTEGRAL_COLUMN_STEP!=4)&&(INTEGRAL_COLUMN_STEP!=2)
#error error defined INTEGRAL_COLUMN_STEP ,must be 2,4,8,16
#endif

#define LOCAL_BLOCK_SIZE ((CL_DEVICE_LOCAL_MEM_SIZE-sizeof(VECTOR_DST)*INTEGRAL_COLUMN_STEP)/sizeof(DST_TYPE)/INTEGRAL_COLUMN_STEP/INTEGRAL_COLUMN_STEP*INTEGRAL_COLUMN_STEP)
#define VECTOR_SRC VECTOR(SRC_TYPE,INTEGRAL_COLUMN_STEP)
#define VECTOR_DST VECTOR(DST_TYPE,INTEGRAL_COLUMN_STEP)
#define VECTOR_INT VECTOR(int,INTEGRAL_COLUMN_STEP)

#define _KERNEL_NAME(s,d,t) prefix_sum_col_and_transpose_##s##_##d##_##t
#define KERNEL_NAME(s,d,t) _KERNEL_NAME(s,d,t)

#if INTEG_TYPE == INTEG_SQUARE
#define compute_src(src) (src*src)
#define _kernel_name_ KERNEL_NAME(SRC_TYPE,DST_TYPE,integ_square)
#elif INTEG_TYPE == INTEG_COUNT
#define compute_src(src) ((DST_TYPE)(0)!=src?(DST_TYPE)(1):(DST_TYPE)(0))
#define _kernel_name_ KERNEL_NAME(SRC_TYPE,DST_TYPE,integ_count)
#elif INTEG_TYPE == INTEG_DEFAULT
#define compute_src(src) src 
#define _kernel_name_ KERNEL_NAME(SRC_TYPE,DST_TYPE,integ_default)
#else
#error unknow INTEG_TYPE by complier with options -D
#endif

#define VLOAD FUN_NAME(vload,INTEGRAL_COLUMN_STEP)
#define VSTORE FUN_NAME(vstore,INTEGRAL_COLUMN_STEP)



inline VECTOR_DST read_vector_from_local(const __local DST_TYPE* local_block,const int column_size,VECTOR_INT* index_v){
	VECTOR_DST v=(VECTOR_DST)(0);
	for(int i=0;i<column_size;++i){
		((DST_TYPE*)(&v))[i]=local_block[((int*)index_v)[i]];
	}
	return v;
}
inline void write_vector_to_local(__local DST_TYPE* local_block,const int column_size,VECTOR_INT* index_v,const VECTOR_DST v){
	for(int i=0;i<column_size;++i){
		local_block[((int*)index_v)[i]]=((DST_TYPE*)(&v))[i];
	}
}
inline void set_vector_elements(VECTOR_DST* v,int from,int to,__global SRC_TYPE* src_ptr){
	DST_TYPE* v_ptr=(DST_TYPE*)(v);
	if(src_ptr)
		for(int i=from;i<to;++i)v_ptr[i]=(DST_TYPE)(src_ptr[i]);
	else
		for(int i=from;i<to;++i)v_ptr[i]=(DST_TYPE)(0);
}

inline void set_block_index(VECTOR_INT* index_v,int index){
	for(int i=0;i<INTEGRAL_COLUMN_STEP;++i,index+=LOCAL_BLOCK_SIZE)((int*)index_v)[i]=index;
}

#define LOAD_VECTOR_SRC(is_aligned,ptr,i)	(is_aligned?( (__global VECTOR_SRC*)ptr)[i]:VLOAD(i,ptr)) 
///////////////////////////////////////////////////////////////////////////////
//! @brief  :	Calculates the prefix sum for each column of an image and transpose
////////////////////////////////////////////////////////////////////////////////
__kernel void _kernel_name_( __global SRC_TYPE *src, __global DST_TYPE * dst, uint width,uint height, uint src_width_step, uint dst_width_step){
 	__local DST_TYPE local_block[INTEGRAL_COLUMN_STEP * LOCAL_BLOCK_SIZE ];	
	__local VECTOR_DST local_block_sum[INTEGRAL_COLUMN_STEP];
	const int local_id=(int)get_local_id(0);
	const int local_size=INTEGRAL_COLUMN_STEP;// the number of local work-items ,equal to get_local_size(0)
	const char V_SHIFT_NUM=7-clz((uchar)INTEGRAL_COLUMN_STEP);	
	const bool is_aligned_src=!(((ulong)src&(INTEGRAL_COLUMN_STEP-1))||(src_width_step&(INTEGRAL_COLUMN_STEP-1)));
	const bool is_aligned_dst=!(((ulong)dst&(INTEGRAL_COLUMN_STEP-1))||(dst_width_step&(INTEGRAL_COLUMN_STEP-1)));
	const int group_step=ALIGN_UP(gf_get_divisor(width,get_num_groups(0)),V_SHIFT_NUM);
	__global SRC_TYPE * src_ptr;
	__global DST_TYPE * dst_ptr;
 	VECTOR_DST last_cycle_sum;
	VECTOR_INT block_index;
	int block_size,column_size,y_step;
	int2 pos;
	for(int group_x=get_group_id(0)*group_step,end_group_x=min((int)width,group_x+group_step);group_x<end_group_x;group_x+=column_size){
		column_size =min((int)(end_group_x-group_x),(int)INTEGRAL_COLUMN_STEP);
		last_cycle_sum=(VECTOR_DST)(0);// reset last_cycle_sum
		for(int y=0; y<height; y+=block_size){
			block_size = min((int)(height-y),(int)LOCAL_BLOCK_SIZE);		
			y_step=gf_get_divisor(block_size,local_size);
			pos=(int2)(group_x,y+y_step*local_id);		
			if(pos.y<height){
				set_block_index(&block_index,pos.y-y);
				VECTOR_DST current_sum = local_id?(VECTOR_DST)(0):last_cycle_sum;
	 			src_ptr	= src + pos.x + pos.y * src_width_step;	 		
				// compute prefix sum in block and save to local memory
				for( int end_index=block_index.s0+min((int)(height-pos.y),(int)y_step);block_index.s0<end_index;++block_index,src_ptr+=src_width_step){
					VECTOR_DST load_src_v;	
					if(column_size<INTEGRAL_COLUMN_STEP){
						if( pos.y+(end_index-block_index.s0)<height-1 ){
							load_src_v=VCONVERT(VECTOR_DST,)(LOAD_VECTOR_SRC(is_aligned_src,src_ptr,0));
						}else{
							load_src_v=(VECTOR_DST)(0);
							set_vector_elements(&load_src_v,0,column_size,src_ptr);
						}
					}else
						load_src_v=VCONVERT(VECTOR_DST,)(LOAD_VECTOR_SRC(is_aligned_src,src_ptr,0));
					current_sum += compute_src(load_src_v);
					write_vector_to_local((__local DST_TYPE*)local_block,column_size,&block_index,current_sum);
				}
				local_block_sum[local_id]=current_sum;
			}			
			barrier(CLK_LOCAL_MEM_FENCE);
			if(!local_id){
				// compute last prefix sum for each block;
				for(int i=1;i<INTEGRAL_COLUMN_STEP;++i)local_block_sum[i] += local_block_sum[i-1];
			}
			barrier(CLK_LOCAL_MEM_FENCE);
			if(pos.y<height && local_id){
				// add last block sum for all current block;
				set_block_index(&block_index,pos.y-y);
				VECTOR_DST tmp;
				for(int end_index=block_index.s0+min((int)(height-pos.y),(int)y_step);block_index.s0<end_index;++block_index){					
					tmp = local_block_sum[local_id-1] + read_vector_from_local((__local DST_TYPE*)local_block,column_size,&block_index);
					write_vector_to_local((__local DST_TYPE*)local_block,column_size,&block_index,tmp);
				}
			}
			barrier(CLK_LOCAL_MEM_FENCE);
			set_block_index(&block_index,block_size-1);
			last_cycle_sum=read_vector_from_local((__local DST_TYPE*)local_block,column_size,&block_index);// save last sum
			// write transposed line to dest matrix		
			if(local_id<column_size){
				dst_ptr	= dst + y + (pos.x+local_id) * dst_width_step; 
				if(is_aligned_dst){
					for(int i=0,end_i=(block_size+INTEGRAL_COLUMN_STEP-1)>>V_SHIFT_NUM,j=local_id*LOCAL_BLOCK_SIZE>>V_SHIFT_NUM;i<end_i;++i,++j)
						((__global VECTOR_DST*)dst_ptr)[i]=((__local VECTOR_DST*)local_block)[j];			
				}else{
					for(int i=0,j=local_id*LOCAL_BLOCK_SIZE;i<block_size;++i,++j)dst_ptr[i]=local_block[j];			
				}			
			}
			barrier(CLK_LOCAL_MEM_FENCE);
		}
	}
}
 